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ABSTRACT 

We test the consistency of active galactic nuclei (AGN) optical flux variability with the 
damped random walk (DRW) model. Our sample consists of 20 multi-quarter Kepler 
AGN light curves including both Type 1 and 2 Seyferts, radio-loud and -quiet AGN, 
quasars, and blazars. Kepler observations of AGN light curves offer a unique insight 
into the variability properties of AGN light curves because of the very rapid (11.6—28.6 
min) and highly uniform rest-frame sampling combined with a photometric precision 
of 1 part in 10® over a period of 3.5 yr. We categorize the light curves of all 20 objects 
based on visual similarities and find that the light curves fall into 5 broad categories. 
We measure the first order structure function of these light curves and model the 
observed light curve with a general broken power-law PSD characterized by a short- 
timescale power-law index 7 and turnover timescale r. We find that less than half the 
objects are consistent with a DRW and observe variability on short timescales (^ 2 
h). The turnover timescale r ranges from ~ 10 — 135 d. Interesting structure function 
features include pronounced dips on rest-frame timescales ranging from 10 — 100 d 
and varying slopes on different timescales. The range of observed short-timescale PSD 
slopes and the presence of dip and varying slope features suggests that the DRW model 
may not be appropriate for all AGN. We conclude that AGN variability is a complex 
phenomenon that requires a more sophisticated statistical treatment. 

Key words: galaxies: active - galaxies: Seyfert - quasars: general - BL Lacertae 
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1 INTRODUCTION 

The continuum of the X-ray through optical radiation spec¬ 
trum observed in active galactic nuclei (AGN) arises in the 
accre tion disk surro u nding the central superm assive black 
hole (IShakural Il973l : IShaGra fc SunvaevI Il973l l. It is well 
known that this region of the spectrum exhibits strong, 
stochastic variability a t the 200 per cent (X-ray) to 10 
per cent (optical) level jMughoj^k y. Done, fc Pounds! [l993l : 
I Wagner fc Witzel 19951 : iKroliklll^^ . Variability is observed 
over a wide range of timescales ranging from hours to months 
to years. The physical mechanisms driving this variability 
are not well understood — m odels range from X-ray flares 
that drive optical variability (iGoosmann et aDl2006ll to lo- 
cal variations in the plasma viscosity of the accretion disk 
llLvubarskiill 199it l caused by small scale angular momentum 
outflo ws triggered by local dynamo processes dKing et al.l 
l2004h . 

The accretion disk of a typical 10 ®Mq supermassive 
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black hole is too small to be imaged. Models of AGN 
variability can only be tested by comparing the stochas¬ 
tic properties of observed AGN light curves to the behav¬ 
ior predicted by the variability model. The most popular 
AGN variability model is the damped random walk (DRW) 
model or 1-parameter Auto-Reg r essive , AR(1), process 
llKellv. Bechtol d. fc Sie m iginow^ 3 |2 009|; iKozlowski et ^ 


I 2 OIOI: iMacLeod 


et al.n2oioin z' 


u et al 


2 OI 3 I I. This model is 


extremely simple and does a good job of fitting ground- 
based variability data, yet we will demonstrate that it fails 
to capture the full range of variability behavior exhibited by 
AGN. We begin by outlining the terminology and mathe¬ 
matics used in discussing variability and the AR(1) process 
in the next section. 


2 THE AR(1) PROCESS 

The light curve of an AGN may be regarded as a discrete 
time series that samples a continuous process. Measurements 
of the source flux Fi are made at times ti for N instants in 
time. Ideally, these measurements are obtained at fixed in- 


© 2015 RAS 



























2 Vishal P. Kasliwal, Michael S. Vogeley, and Gordon T. Richards 


crements separated by a constant sampling interval 5t. In 
such cases, the time interval At between any two flux mea¬ 
surements, Fi & Fj, is related to the sampling interval by 

At = n5t, (1) 


where n = i — j with 1 ^ n < N. This kind of sampling 
pattern is usnally not possible in the case of ground-based 
astronomical time series because of interruptions caused by 
factors such as the weather and the availability of the target 
in the night sky. Even in the case of space-based measure¬ 
ments of the series, it is possible to have ‘missing values’, 
i.e. values of ti with i = nSt, for which there is no corre¬ 
sponding measurement of the flux Fi. The presence of these 
missing values complicates the analysis process. In the case 
of a time series that exhibits stochastic behavior, the goal 
of time series analysis is to characterize the joint probability 
distribution of Fi for the measured data points. If the joint 
probability distribution of the data points is independent of 
time, the time series is said to be stationary, i.e. one would 
observe the same light curve (from a statistical standpoint) 
at all times. 

The timescale on which the bulk fueling rates of the 
AGN vary is in the hundreds of thousands to millions of 
years and therefore cannot be directly probed by obser¬ 
vation for any individual AGN. It is reasonable to expect 
that the stochastic variability observed in the flux time 
series of individual AGN is caused by frequently reoccur¬ 
ring short-lived processes local to the accretion disk, i.e. 
the time series is stationary on the timescales we probed. 
It is well known from the theory of time series analysis 
that a stationary time series can be modeled by an Auto- 
Regressive Moving Average, or ARMA, process jHamiltonl 
1994 Woodward. Gray. fc Elliod 2OI2I: Prado k, We^ gOld : 


Box. Jenkins, fc Reinsell l2006l : iBrockwell &: DavisI l2002l l of 
appropriate order (p,q). The general form of an ARMA(p,q) 
process, i.e. a process with p autoregressive terms and q 
moving average terms, is given by 


Fi — ^ ( (j^mFi — jn T ^ ^ ^i — nnJi — n -f Wi, 


( 2 ) 


where (pm are autoregression coefficients. On are moving av¬ 
erage coefficients, and Wi are known as ‘innovations’ in sta¬ 
tistical literature. Autoregressive terms introduce infinite- 
duration correlations in the data while moving average terms 
introduce finite-duration correlations. It is possible to repre¬ 
sent any stationary process as a pure autoregressive or mov¬ 
ing a verage process, albeit perhaps of infinite order (IScargld 
[1983). Innovations drive the stochastic behavior of the pro¬ 
cess. Each innovation is an independent random number typ¬ 
ically drawn from a Gaussian random distribution with zero 
mean and variance P oc St. To make the process station¬ 
ary we must restrict \(j)i -ni\ < 1 for all values of m. It is 
clear from previous work (iK ellv. Bechtol d. fc Siemiginowskal 
l2009l : lKozlowski et al.l[^010l : lMacLeod et al.ll201ol l that AGN 
light curves exhibit long-term correlations. The simplest 
ARMA process that exhibits long-term correlations is the 
purely autoregressive (no moving average terms) AR(1) pro - 
cess or damped random walk (DRW) dlvezic et al.l 120141 ') 
given by 


Fi = (piFi-i -I- Wi, 


(3) 


with 0 < 01 < 1. The lone autoregression coefficient 0i con¬ 


nects the current value of the series to the previous value, 
making the AR(1) process a memoryless or Markov process. 
The numerical value of 0i depends upon the sampling in¬ 
terval St. For a given sampling interval, a value of 0i ~ 1 
implies that Fi ~ Fi-i and so the value of 0i sets a char¬ 
acteristic timescale r after which the contribution of the 
previous value of the series will be relatively insignificant, 
i.e. T is the decorrelation timescale. Numerically, 

_ 

01 = e - , (4) 

where St is the sampling interval and r is the aforementio ned 
characteristic decorrelation timescale dGillespiel Il996h . If 
St 2^ T, then 0 ~ 0 and any given value of the time se¬ 
ries is only very weakly correlated with the previous value. 
On the other hand, if St t, then 0 ~ 1 and the time series 
exhibits strong correlation. Eor the series to be stationary, 
we must have 0i < 1, i.e. St < t. It turns out that the vari¬ 
ance of the innovations is also linked to the decorrelation 
timescale: 

(wi) (5) 

where a allows the amplitude of the variance of the innova¬ 
tions to vary independently of St and r. Eor a fixed value of 
a, if the sampling interval is much longer than r, the vari¬ 
ance of the innovations will be very close to a. Combined 
with the accompanying small value of 0i, the resulting time 
series will look close to pure white noise. If, on the other 
hand, we choose a sampling interval much smaller than the 
decorrelation timescale, the variance of the innovations is 
forced to be very close to zero and the value of 0i is very 
close to 1, resulting in very obvious correlations between the 
data points. 

Physically, the AR(1) process could result from an ac¬ 
cretion disk with local ‘spots’ that contribute more or less 
flux than the mean flux level of the disk. These spots ap¬ 
pear at rando m and dissipate over s ome characteristic phys¬ 
ical timescale (lAgol fc Dexterll201lll . Under the assumption 
that the distribution of the contribution to the total flux 
from each spot is Gaussian, the Central Limit Theorem as¬ 
sures us that the overall changes in flux, i.e. the innovations 
in Eq. m and in Eq. ®, are drawn from a Gaussian dis¬ 
tribution. Long term correlations exist because the spots 
do not dissipate instantaneously. The power spectral density 
(PSD) of this time series follows a bent power-law given by 

while the auto-covariance function (ACVF), related to the 
PSD via a Fourier Transform (Wiener-Khintchine theorem), 
is given by 

ACVFAR(i){^t)=Pe-^, (7) 

where a is the low-frequency asymptotic amplitude and r is 
the cha racteristic time scale at which the PSD turns over and 
decays (lGillesDielll996ll . Hence the auto-correlation function 
(AGF) is given by 

ACFAR(i){M) = e-^, ( 8 ) 

Does the AR(1) process actually descri be the variabil¬ 
ity of AGN? Recen t ground based studies dZu et al.ll2013l : 
iGraham et al.ll2014l l indicate that on very short timescales. 
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AGN light curves show strong autocorrelation with PSD 
slopes steeper than 1//^. To model stochastic light curves 
with arbitrary PSD slopes, we introduce the damped power- 
law (DPL) model by generalizing the AR(1) PSD in Eq. ([ 6 |) 
to 


PSDdplU) 


„2 

'^Eff 

1 + (27r/T)T'’ 


(9) 


where we have lumped the quantities in the numerator of 
Eq. ® into a single variable ctb// and introduced 7 to al¬ 
low the logarithmic slope of the PSD to be a free parameter. 
When 7 < 2, the process exhibits weaker autocorrelation on 
short timescales than the AR(1) i.e. the time series looks 
less smooth. When 7 > 2, the process shows stronger au¬ 
tocorrelation on short timescales than the AR(1) i.e. the 
time series looks smoother. This DPL model is a simplifi- 
cation of the popular ‘Bending Power-Law’ (BPL) model of 
iMcHardv et al.l ll2004h : 


L 

PSDBPL{f) = Af-‘^^l[ 

i=l 



( 10 ) 


where Oi are the logarithmic PSD slopes above the corre¬ 
sponding ‘bend’ frequencies fsi- Allowing the logarithmic 
PSD slope to be a free parameter is similar to the s tructu re 
function parametrization chosen bv ISchmidt et all (l201Cll 'l: 


SF{At) = A 



( 11 ) 


where A plays the same role as 7 in the DPL model-they 
allow the correlation level at fixed time-lag to vary from that 
predicted by a pure AR(1) process. To test our DPL model, 
we use light curve data from the NASA Kepler mission. 


3 THE KEPLER MISSION 

The Kepler mission was designed to survey the Miky Way 
for Earth-sized and smaller planets around the habitable 
zone by using the transit-method to detect exo-planets. 
Kepler is designed as a Schmidt-camera with a primary 
aperture of 0.95-m, a fixed 115.6 deg^ FOV, and a half¬ 
maximum bandpass from 435-nm to 845-nm with ^ 1 
per cent relative spectral response from 420 nm to 905 
nm dKepler Instrunient Handbook! [kSCI- 19033I1 . The Ke¬ 
pler focal-plane consists of 42 CCDs arranged in 21 mod¬ 
ules that are readout by 4 channels each. Each 50 x 25-mm 
CCD has 2200 x 1024 pixels yielding an image scale of 3.98 
arcsec per pixel. The Kepler PSF is intentionally large-6.4 
pixels near the center of the FOV-in order to increase the 
photometric precision dCillilandl l2004l ~) by ensuring that no 
pixel contains more than 60 per cent of the flux from a 
point source target. This has two consequences that com¬ 
plicate simple aperture photometry: 1 . point sources, such 
as stars, appear extended, and 2 . faint background sources 
‘crowd’ the aperture used to measure the flux from a tar¬ 
get of interest. The signal at each pixel on the CCDs is 
measured over an interval known as the integration time or 
frame which consists of a variable exposure time and a fixed 
readout time of 0.520 s. The flight exposure time is 6.02 s 
yielding a flight integration time of 6.54 s per frame. Two 
datasets can be created by Kepler-, the long eadence (LC) 
data and the short cadence (SC) data. The long cadence 


data is integrated over 270 frames (29.430 min), while the 
short cadence data is integrated over 9 frames (0.981 min). 
Most of the objects observed by us have no SC data-only LC 
data was obtained. Furthermore, the SC data set is harder 
to calibrate b ecause of the small nu mber of SC quiet tar¬ 
gets observed llKinemuchi et al.ll2012ll. For these reasons we 
only use LC data in our analysis. Kinemuchi et al.l il2012l ~l 
provides an excellent overview of the final data products 
supplied by Kepler, including strategies for dealing with the 
various types of artifacts present in the data. 

The LC data-set poses several calibration chal¬ 
lenges. Certain readout channels on specific CCD mod- 
ules have known performance issues di scussed in the 
iKepler Instrument Handbook! ||KSCI-19033I1 including out- 
of-spec read noise levels and gain. More troubling are the 
‘ moire’ and ‘rolling band ’ effects seen in some channels 
llKolodzieiczak et al.11201^ '). The combination of module and 
channel is usually quoted in the format module^^^:.channel#. 
A Kepler target will fall on a fixed sequence of 4 module and 
channel combinations over the course of a year, i.e. since the 
spacecraft rolls 90° every quarter, at the beginning of every 
5*^ quarter targets will fall on the same module#.channel# 
combination and will be subject to the same instrumenta¬ 
tion issues. We discuss the impact of these effects in Sec. 

El 


Spacecraft operation events and systematic trends 
such as thermally-driven focus variations, pointing offsets, 
and differential velocity aberration (DVA), introduce 
artifacts into the data; further systematics are intro¬ 
duced as the result of light losses caused by the time 
dependent sampling of the outer wings of the tar- 
get PSF and contamination from n ei ghboring sources 
aKepler Data Characte ristics Handbook] | kSCI- 19040-()03 : 
iKepler Instrument HandbookI IkSCI- 190331 . Such artifacts 
can be misinterpreted as being astrophysical in origin, and 
can mask true astrophysical signals. Post-calibration Kepler 
data is available in both uncorrected (SAP flux) and cor¬ 
rected (PDCSAP flux) forms. Corrections to the calibrated 
data are obtained usin g the PDC module (|Stumpe et al.l 
l 2012 l:ISmith et al.ll 2012 ll . 

The PDC module identifies features common to hun¬ 
dreds of quiet targets on each CCD by examining the cross¬ 
correlation between targets. The most correlated targets are 
used to create 16 best-fit vectors called Cotrending Basis 
Vectors (CBVs) using Singular Value Decomposition (SVD). 
Corrections to the calibrated data are obtained by remov¬ 
ing the CBVs from the data using a weighted normaliza¬ 
tion. The algorithm used by PDC to compute these CBV 
weights changed with the introduction of Ver 8.0 of the Ke¬ 
pler Science Operations Center (SOC) pipeline (Data Re¬ 
lease 14). Prior to pipeline Ver. 8.0, the PDC module used 
a Least Squares (PDC-LS) approach to compute the CBV 
weights. Ver. 8.0 introduced a Bayesian Maximum A Poste¬ 
riori (PDC-MAP) algorithm compute the CBV weights and 
optimally remove spacecraft-induced trends while minimiz¬ 
ing t he removal of the true underlying signal dStumpe et al.l 
I2OI2II . The Bayesian-MAP procedure assumes that the the 
observed signal may be represented as 


F = He-\-e, 


( 12 ) 


where 77 is a design matrix consisting of the CBV vectors 
(typically 4), 6 is the vector of the CBV weights (which 
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is what the PDC module is trying to determine), and e ~ 
(0, Ejv) is the (uncorrelated, i.i.d.) instrument noise vector. 
The log-likelihood of an observed light curve, assuming this 
model for the data, is then given by 

lnp{F\e) = + ln(|S^|) + {F - - HO) 

^ (13) 

Maximization of this likelihood can be performed analyti¬ 
cally, however the resulting de-trended light curve will nec¬ 
essarily have intrinsic variability removed because LS cal¬ 
culates the CBV weights by projecting the raw light curve 
onto each CBV—any coincidences between the true light 
curve and the CBV will spuriously add to the CBV weight. 

To avoid this behavior, the new PDC-MAP approach 
applies a weighted Bayesian prior, p{6), to the likelihood 
in to bias the best-fit CBV weights toward removing 
only the spacecraft-induced variability signal while leaving 
intrinsic variability intact. p{0), is constructed by using LS 
to de-trend quiet targets in the phase-space vicinity (RA, 
Dec, and apparent magnitude) of the object of interest. The 
prior weight, Wpr, is calculated such that it is closer to 0 for 
relatively quieter targets and closer to 1 for relatively active 
targets. The PDC algorithm computes 

9 = argmax(lnp(T’|0) -|- Wprlnp(0)). (14) 

e 

Therefore p{9) biases the CBV weights towards values cor¬ 
responding to those for neighboring quiet targets. Since the 
CBV weights are biased towards removing just enough fea¬ 
tures in the light curve to render featureless light curves for 
neighboring quiet targets, intrinsic variability in the light 
curve of a truly variable object should be relatively unaf¬ 
fected by the de-trending procedure. Further d etails on the 
PDC de-trending algorithm may be found in ISmith et al.l 
1120121 ) . 

The benefit of using the Aep/er-MAST supplied PD- 
CSAP flux light curves is that for most targets, the PDC- 
MAP algorithm produces optimally de-trended light curves 
free of instrument effects. At the same time, the drawback 
of doing so is that it is theoretically possible for the PDC- 
MAP pipeline to remove intrinsic variability, particularly 
when the target exhibits intrinsic long-term trends that can¬ 
cel out the systematic trend, or when the location of the 
target in position-luminosity phase-space is lacking an ade¬ 
quate number of close neighbors. On the other hand, using 
the raw SAP flux light curves and performing a customized 
de-trending using the PyKE tool kepcotrend implies that 
the CBV weights must be computed using no information 
from neighboring targets and hence are likely to project out 
intrinsic variability in the light curve. 

We use Kepler-MAST supplied light-curves from Data 
Release 23 that have been processed using version 9.1 of the 
Kepler SOC pipeline (PDC module uses PDC-MAP algo¬ 
rithm) and only perform a minor correction to remove the 
large quarterly offsets in the data caused by spacecraft rolls. 

4 SAMPLE SELECTION 

Only ~ 7 AGN were known to lie in the Kepler FOV prior to 
the start of the mission due to the lack of coverage from ex¬ 
isting extragalactic surveys. Since the beginning of the mis¬ 
sion in May 2009, Guest Observer (GO) efforts to find more 


AGN in the Kepler FOV by multiple groups dCarini fc Rvld 
I2OI2I : lEdelson fc Malk^l2012l : IWehrle et al.ll2013h have led 
to a total of roughly 80 AGN of various types being mon¬ 
itored by Kepler. The May 2013 failure of a critical reac¬ 
tion wheel led to the termination of the original Kepler mis¬ 
sion. We wish to examine variability properties across a wide 
range of AGN type and redshift; however a large number of 
the Kepler AGN were selected photometrically and do not 
have redshifts. As such, we restrict our analysis to the 20 
that are spectroscopically confirmed AGN. This allows us to 
reject non-AGN contaminants, and also allows us to study 
variability behavior as a function of AGN type and perform 
comparisons in the rest-frame of the object. 

Of the 7 AGN known to lie in the Kepler FOV prior 
to 2009, the best studied is kplr006932990, also known as 
Zw 229-15. The VGV Catalog llVeron-Cettv fc Veronllioloh 
lists kplr006932990 as a radio-quiet Type 1 Seyfert at a red- 
shift of 0.027 with V-band absolute magnitude My = —19.9. 
Based on H/3 reverberation mapping. iBarth et al.l (1201 ll) ob¬ 
tain a virial black hole mass of Mbh = 1.00j)Q'24 x IO^Mq. 
The other Type 1 AGN known to exist in the Kepler FOV 
prior to the commencement of the mission is the X-Ray 
source kplr09650715 (RXS J19298-I-4622), which is a radio¬ 
quiet Seyfert 1 at z = 0.127 with V-band absolut e mag¬ 
nitude Mv = —21.8 (IVeron-Cettv fc VeronI I2OI0I) . Two 
narrow-line objects were known to exist in the Kepler FOV: 
kplr00870 3536 (IGR J19 473-k4452), which is a radio-quiet 
Seyfert 2 (iMasetti et ahlliood) at z = 0.0539 with V-band 
absol ute magnitude Mv = —21.1 dVeron-Gettv fc VeronI 
I2OI0I) . and kplr011021406 (6G 1908-1-4829), which is a radio- 
loud object at 2; = 0.513 with photographic magnitude 
Mq = —21.9 variously class ified as either a S eyfert 1.5 by 
IVeron-Gettv fc VerorJ ll2010l ) or a FSRQ by iHealev et al.l 
1 2OO7I) . At higher redshifts. pre-Kepler AGN include the 
FSRQ blazars kplrll606854 (MG4 J191843-I-4937) at 2 = 
0.926 and k p lr010 663134 (Q 1922+4748) at 2 = 1.52 
iHealev et al.l (l2007l ). The last of the pre-Kepler AGN is 
the quasar kplr012208602 (4G 50.47) at a redshift of 1.098. 
Spectra obtained using the MMT indicates that this ob¬ 
ject has an exceptionally strong GIV 1549 line with a rest- 
frame equivalent width of 240 A and a power-law spectrum 
from 178 MHz to 5 GHz with spectral index a = —0.65 
dWalsh et al.lll983). Both kplr01 2208602 and kplr010663134 
possess radio-jets dLiu fc Zhanglliool ). 

lEdelson fc Malk^ d2012l) obtained spectra of the re¬ 
maining 13 AGN in our sample usin g the Kast double 
spectr ograph on the Lick 3 m telescope. lEdelson fc Malkanl 
d2012l) classify kplr006690887 as a BL Lac object based 
on a featureless spectrum and coincidence with the radio 
source NVSS J192631+420958, and classify kplr09215110 
as a Type 1.9 Seyfert. Since no sub-classification is pro¬ 
vided for the remaining 11 objects, we performed a visual 
inspec tion of the spectra provided by lEdelson fc Malkanl 
d2012h to characterize these objects by emission line type. 
We conclude that while most of the remaining objects 
are Type 1 AGN with prominent broad emission lines, 
kplr003347632 may be a Type 2 object. Although no SDSS 
spectra exist for the AGN in our sample, 4 of the new 
AGN (kplr002694186, kplr002837332, kplr003337670, and 
kplr005597763) land within the SDSS DRIO footprint. By 
combining SDSS with STScI DSS imaging, we sub-classify 
the objects as Seyferts or QSOs based on visual appearance 
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Table 1. Kepler AGN Sample. 


KeplerlD 

Identifier 

Alt. ID 

RA 

Hrs min sec 

Physical Parameters 

Dec z Mg 

Deg min sec 

AGN Type 

6932990 

Zw229-15l' 

19 05 25.969 

-1-42 27 40.07 

0.0275 

-23.8 

sJTf 

12158940 

1925-I-50* 

19 25 02.181 

-1-50 43 13.95 

0.067 

-21.9 

Syl 

11178007 

W2R 1858-I-48* 

18 58 01.111 

-1-48 50 23.39 

0.079 

-20.2 

Syl 

9650715 

RXS J19298-I-4622T 

19 29 50.490 

-1-46 22 23.59 

0.127 

-21.8^ 


2837332 

W2R 1910-I-38* 

19 10 02.496 

-1-38 00 09.47 

0.130 

-20.4 

Syl 

3347632 

W2R 1931-I-38* 

19 31 15.485 

-1-38 28 17.29 

0.158 

-21.0 

Sy2? 

5781475 

W2R 1915-I-41* 

19 15 09.127 

-1-41 02 39.08 

0.222 

-22.1 

Syl 

2694186 

W2R 1904-I-37* 

19 04 58.674 

-1-37 55 41.09 

0.089 

-23.8 

SyT’ 

9215110 

W2 1922-I-45* 

19 22 11.234 

-1-45 38 06.16 

0.115 

-22.2 

Syl.9* 

10841941 

W2R 1845-I-48* 

18 45 59.577 

-1-48 16 47.57 

0.152 

-21.2 

Syl- 

3337670 

W2R 1920-I-38* 

19 20 47.750 

-1-38 26 41.28 

0.368 

-23.0 

Syl 

7610713 

1931-I-43* 

19 31 12.566 

-1-43 13 27.62 

0.439 

-24.7 

QSO 

6690887 

W2R 1926-I-42* 

19 26 31.089 

-1-42 09 59.12 

0.154 

-21.7 

BL-LaC’!* 

6595745 

W2R 1914-I-42* 

19 14 15.492 

-1-42 04 59.88 

0.502 

-24.6 

QSO 

5597763 

W2R 1853-I-40* 

18 53 19.284 

-1-40 53 36.42 

0.625 

-25.1 

QSO 

11606854 

MG4 J191843-I-4937’ 

19 18 45.617 

-1-49 37 55.06 

0.926 

-25.0 

FSRQt 

10663134 

Q 1922-|-4748t 

19 23 27.234 

-1-47 54 17.00 

1.520 

-25.2 QSOt; FSRQt; JetH 


8703536 IGR 319473+4452'^ 19 47 19.308 +44 49 42.32 0.0539 -21.7 

11021406 6C 1908+4829^ 19 09 46.501 -h48 34 32.26 0.513 -23.2 Syl.5‘1'; FSRQ^ 

12208602 _ 4C 50.47t 19 26 06.318 -h50 52 57.14 1.098 -24.7 QSPt; JetH _ _ 

* Reliable Identifications of Active Galactic Nuclei from the WISE^ 2MASS, a nd ROSAT All-Sky Surveys dEdelson MalkarjF2012t) 
^A Catalog of Quasars and Active Nuclei: 13th Ed. llVeron-C etty &: Veron|| 2010ll 

^ An All-Sky Survey of Flat-Spectrum Radio Sources l|Healev e^alH^OQTll _ 

^Kepler Observations of Rapid Optical Variability in the BL Lac Object W2R1 926H-42 Isdelson et al.|[2013h 

IIA New List of Extra-Galactic Radio Jets l|Liu &; Zhandl20Q2ll _ 

^CGRaBS: An All-Sky Survey of Gamma-Ray Blazar Candidates IIHealev et al.ll2008h 
’Weak spectral features-unsure sub-classification 

The AGN are grouped based on the light curve categories discu ssed in Sec. Col. 6 l i sts SP SS g-band absolute magnitudes computed 
from the SDSS apparent magnitude supplied in the lKepler Input Catal^ ll201ll) with cosmologic al parameters 
Hq = 67.77 kms~^Mpc~^, SJjvf = 0.3071, and Qa = 0.6914. We performed a k-correction using Eq. (10.29) in IPetersonI ll2003li assuming 
OL = —0.5 llRichards et al.ll2006h to compute the absolute magnitude at 2 = 0. W e present tentative s ub-cla ssifications for 11 AGN in 
Col. 7 using imaging from the STScI Digitized Sky Survey, SDSS and spectra in lEdelson fc MalkanI <l2012t) . In the case of broad-line 
objects we distinguish between Seyfert Is and QSOs based on the presence of a galactic component in STscI DSS/SDSS imaging and 

also based on Mg > — 24 for Seyferts. 


and k-corrected ab solute magnitude (tra ditionally Mi < 
—22 for QSQs from iR.icharc s et ^ (j2006l )) computed from 
the iKepler Input Catalog! ( ' 2 OIII) SDSS g-band apparent 
magnitude available on MAST. 

Table m presents positional and physical parameters 
for the light curves from Fig. [S] Column 1 lists the Ke- 
plerlD used when querying the Kepler M AST database 
(jKepler Archive Manuall KDMC- 10008-004 ) . Column 2 lists 
the ‘common’ name of the object along with a pri- 
mary reference for the obje ct, usually drawn from eithe r 
IVeron-Cettv fc VeronI (l201(]l l or lEdelson fc MalkanI ll2012ll . 
The next two columns supply the J2 000.0 RA and Dec lo¬ 
cation of the object as listed in the iKenler Input Catalod 
(120111) . The next column supplies the redshift of the ob¬ 
ject as listed in the primary reference for the object. The 
next column supplies the k-corrected SDSS g-band absolute 
magnitude [Mg) based on the appar ent SDSS g-band KIC 
magni tude {nig). As described in the lKepler Input Cataloa 
(I2OIII) , the galactic extinction corrected rrig was determined 
during the Stellar Classif ication Project (SCP) by using 
DAOPHOT (IStetsonlll987l ) to perform PSF photometry on 
star-like objects. We have converted these g-band magni¬ 
tudes to_absolute_magnitudes at z = 0 using k-corrections 
from IPetersoiil (l2003l ) with cosmological parameters Hq = 


67.77 kms“^Mpc“^, = 0.3071, a nd Da = 0.6914, and 
assuming a = —0.5 (iRichards et al.ll2006l ). Column 7 lists 
the AGN type of the object either from the primary refer¬ 
ence or based on a visual exami nation of imaging and spec¬ 
tra from publicly archive d data (iKepler Input Catalodl^lll : 
lEdelson &: Malkanll2012l ). 

Light curves from Kepler exhibit huge quarterly off¬ 
sets and require ‘stitching’ before we may analyze them. 
The next section discusses how we stitch the light curves 
together. 


5 STITCHING LIGHT CURVES ACROSS 
QUARTERS 

As discussed in Sec. [3] and U] quarterly discontinuities ex¬ 
ist in Kepler data. Once every 3 months, the spacecraft 
performs a roll maneuver to transmit data to the Earth. 
After each roll maneuver, targets fall on different CCDs re¬ 
sulting in slightly different distributions of the target flux 
across adjacent pixels. To compensate, the target aperture 
is redefined after every roll maneuver. Due to the large size 
of the Kepler PSF relative to the aperture size, different 
apertures contain different portions of the flux from the tar- 
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get. This results in the severe discontinuities observed in 
the target flux. Fig. [T] shows the discontinuities present in 
the PDCSAP light curve of kplr006932990. The top panel 
is plotted in the observed frame where the sampling inter¬ 
val is 29.42 min. The bottom panel is plotted in the rest 
frame of kplr006932990 with sampling interval 28.64 min 
because of the cosmological redshift factor of 1 -I- z = 1.0275 
in the case of this object. The effect of this cosmological 
time-contraction is to reduce the effective Kepler sampling 
interval for more distant objects resulting in shorter but 
more densely sampled light curves. The discontinuities in 
the top panel of Fig. [T] must be removed before the long 
term variability properties of this object can be studied. One 
could remove such discontinuities by re-calibrating the data 
on an object-by-object basis, but this requires considerable 
manual effort. Kepler data products include both the light 
curves of the target as well as the Target Pixel File (TPF) 
corresponding to the target. The TPF consist of ‘postage 
stamp’ snapshots of the object taken at every cadence. The 
light curve of the object is constructed by defining an aper¬ 
ture consisting of a subset of the pixels in the postage 
stamp. This aperture is determined by th e Kepler pipeline 
based on the lKepler Input Catalog (|201l[) and is optimized 
for stellar targets as outlined in [Data Processing HandbookI 
(IKSCI-19081-001|) . By re-defining the target aperture by 
hand, it is possible to reduc e the discontinuities in the flux 
across quarters as shown in ICarini fc Rvlj ll2012ll . Such re¬ 
extracted light curves must then be de-trended of spacecraft 
induced features by removing the CBVs (discussed in[3| with 
custom weights determined using the PyKE tool kepcotrend. 

A simpler but cruder approach to removing quarterly 
discontinuities is to match a suitable metric across quar¬ 
ter bo undaries dKinemuchi et all l2012l l. iMushotzkv et all 
(l201llf perform a simple end-matching to stitch light curves 
across quarters, i.e. they correct the measured flux values in 
the second of every sequential pair of quarters by the differ¬ 
ence between the flux value of the last point of the leading 
member, and the flux value of the first point of the trail¬ 
ing member of the pair of quarters. This method is quite 
su itable for the high sign al-to-noise (S/N) objects studied 
in iMushotzkv et al.l (l201ll '). but does not work as well as the 
variability signal level drop s closer to the i nstru ment noise 
level in the dimmer targets. lEdelson et al.l (l2013l l stitch the 
light curve of the BL Lac object kplr006690887 across quar¬ 
ters by using the per-quarter average flux value to determine 
a multiplicative facto r that they then use to scale quarters. 
iRevalski et al.l (120141 ') use a similar procedure but determine 
their multiplicative scaling factor from the average of the 
last 20 points in the quarter preceding the discontinuity and 
the first 20 points in the quarter following the discontinuity. 

We stitch the light curves across quarterly discontinu¬ 
ities by matching the average flux value of the last 100 points 
of the leading quarter to the first 100 points of the trailing 
quarter in every sequential pair of quarters. The method is 
robust to outliers and is relatively unaffected by low S/N. It 
may be argued that it is more correct to use a fixed duration 
in time at the beginning and ends of quarters to compute 
average fluxes, rather than a fixed number of points. Using 
a variable number of points determined by a constant time 
window introduces variations in the reliability of the statistic 
for removing discontinuities; a one day window in the refer- 




Figure 1. PDCSAP light curve of the Seyfert 1 kplr006932990 
(Zw 229-15) — Top panel: unstitched; bottom panel: stitched. 
The huge discontinuities in flux in the top panel are caused by 
redefinition of the target aperture after the spacecraft performs 
a roll maneuver as explained in Sec. [5] and are endemic to all the 
Kepler light curves. Our stitching technique serves to remove the 
quarterly discontinuities present in the light curve obtained from 
MAST. 


ence frame of a a = 1 object will have twice as many flux 
measurements than the same window in the case of a 2 = 0 
object. Fixing the number of points has the benefit that it 
is equally applicable to stitching together the light curves 
of objects at unknown redshifts - the stitched light curves 
of such objects can stand on equal footing with the stitched 
light curves of objects with known redshift. Fig. [T] shows an 
application of our method to the high S/N light curve of 
kplr006932990. The top panel presents the un-stitched light 
curve while the bottom panel presents the stitched light 
curve corrected to the rest-frame of the AGN. As can be 
seen in the figure, our stitching method removes the dis¬ 
continuities present in the original light curve. Similar re- 
sults may be exp e cted from the end-matching technique of 
IMushotzkv et all (1201 ll ). In contrast to the high S/N light 
curve of kplr006932990. Fig. [2] shows the un-stitched and 
stitched light curves of the low S/N AGN kplr003337670. 
The relative thickness of this light curve (caused by scat¬ 
ter from measurement noise) makes it essential to consider 
a suitably large number of points at the ends of quarters 
when stitching. Simple end-matching fails to stitch this light 
curve, and using a smaller number of points in the stitching 
results in obvious flux mis-matches. We use our stitching 
algorithm to remove the quarterly discontinuities in all 20 
AGN light curves and discuss the resulting multi-quarter 
Kepler light curves in the next section. 
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Figure 2. PDCSAP light curve of the Seyfert 1 kplr003337670 — 
Top panel: unstitched; bottom panel: stitched. The noisiness of 
this light curve makes it essential to use a large (100) number of 
points at the edges of each quarters when performing the stitch. 
End-matching is completely ineffective while using a smaller num¬ 
ber of points at the quarter edges results in obvious mis-matches 
across the quarter boundaries. 


6 AGN LIGHT CURVE FEATURES 

We divide the AGN in the sample into 5 categories of ob¬ 
jects based on the visual appearance of their light curve. Fig. 
[3] shows all 20 AGN grouped by visual features present in 
each light curve. Within each category, we order the AGN 
by 2 . All 20 light curves are plotted over the same rest-frame 
time range on the x-axis. Category 1 consists of 3 objects 
that have the most stochastic-looking light curves i.e. light 
curves bereft of any sinusoidal features or flares. Weak rip¬ 
ple features appear in the light curves of Category 2 objects 
(4 AGN). Category 3 consists of 5 objects that have pro¬ 
nounced ripple features in their light curves, while category 
4 objects have light curves with flaring behavior (5 AGN). 
The last category of 3 objects have primarily featureless light 
curves consistent with marginal levels of variability. 

The light curves show a wide variety of different types of 
behavior. There are indications that individual light curves 
can show more than one type of feature. Category 1 ob¬ 
jects (kplr6932990 , kplr012158940, and kplr011178007) are 
shown in the U* row of Fig. |3] These objects have the most 
‘stochastic’-looking light curves and are mostly free of any 
recurring features and trends. They are also the closest AGN 
in this study (2 < 0.1), are all Seyfert Is, and have the 
largest amplitudes in variability amongst all of the AGN in 
this study barring only kplr006690887. 

Four AGN (kplr009650715, kplr002837332, 

kplr003347632, and kplr005781475) have mild oscilla¬ 
tory features in otherwise stochastic-looking lightcurves. 
This behavior is a mixture of the behaviors exhibited by the 


Table 2. Kepler AGN Detector Sequences. 


Identifier 
Kepler ID 
6932990 
12158940 
11178007 
9650715 
2837332 
3347632 
5781475 
2694186 
9215110 
10841941 
3337670 
7610713 
6690887 
6595745 
5597763 
11606854 
10663134 
8703536 
11021406 
12208602 


Detector Properties 
Module#. Channel# 
14.48*, 8.24, 12.40*, 18.64 
24.81*, 10.29, 2.1, 16.53 
23.79*, 15.51§, 3.7’, 11.35§ 
14.46, 8.22t, 12.38, 18.62* 
22.74, 20.70, 4.10ll, 6.14* 
16.55, 24.83*, 10.31, 2.3 
7.17, 17.57, 19.65, 9.25 
22.75*, 20.71, 4.11*, 6.15 
13.41, 13.42, 13.43, 13.44** 
20.71, 4.11*, 6.15, 22.75* 
16.53, 24.81*, 10.29, 2.1 
8.22*, 12.38, 18.62*, 14.46 
12.37, 18.61, 14.45, 8.21 
8.22*, 12.38, 18.62*, 14.46 
23.80, 15.52*, 3.8’, 11.36* 
24.82, 10.30*, 2.2, 16.54 
19.66*, 9.26*, 7.18, 17.58** 
9.27, 7.19, 17.59, 19.67 
23.77, 15.49, 3.5’, 11.33 
24.81*, 10.29, 2.1, 16.53 


*‘medium’ Rolling Band & Moire 
'*^‘high’ Rolling Band Sz Moire 
^Out of Spec Read Noise 
§Out of Spec Undershoot 
II Out of Spec Gain 
^CCD Failed - no data 


The AGN are grouped based on light curve features as in Ta- 
ble^ Col. 2 lists the module#.channel# sequence for each AGN. 
Note that some objects share this sequence of module#.channel# 
though the starting values may be different. This implies that 
such AGN are affected by similar levels of instrumentation arti¬ 
facts and is used in Sec. |6]to decide that instrumentation is not 
responsible for some of the features observed in the light curves. 


category 1 AGN above and the category 3 AGN discussed 
below. The light curves of these AGN are shown in the 
2"“^ row of Fig. |3] between the light curves of the category 
1 and 3 AGN to facilitate comparisons. Although there 
are indications of oscillatory behavior in this intermediate 
category of objects, the oscillations have very poorly defined 
time-periods as compared to objects in category 3. However, 
their light curves are not as purely stochastic looking as the 
light curves of the objects in category 1. Supporting this 
notion of an intermediate state is the observation that the 
object kplr012158940, which is grouped with the category 
1 AGN, appears to be in the process of switching between 
categories as it begins to display features more reminiscent 
of category 3 light curves toward the end of the data. 

Category 3 AGN (kplr002694186, kplr009215110, 
kplr010841941, kplr003337670, and kplr007610713) are 
shown in the 3'^'^ row of Fig. [3] These objects appear to 
exhibit pronounced rippling features in the light curve. 
The aforementioned ‘moire’ and ‘rolling-band’ effects that 
are known to exist in specific detector modules and chan¬ 
nels in Kepler are the natural suspects. These effects are 
known to occur on timescales of a few days and have pat- 
terns very similar to thos e observed in these light curves 
llKolodzieiczak et al.l 120101) . However, a close examination 
of the PDCSAP data suggests that the oscillations observed 
in these light curves are seen even when the target lands on 
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Figure 3. ‘Stitched’ light curves of all the AGN in our investigation. The AGN are grouped based on the light curve categories discussed 
in Sec.|6] Light curves within a category are sorted in order of increasing redshift. All 20 light curves are plotted over the same rest-frame 
time range on the x-axis. Refer to Table [J for details on individual objects. 


a detector module and channel combination that is known 
to be free of the moire and rolling-band effects. For exam¬ 
ple, kplr002694186 exhibits strong oscillatory features in the 
light curve as seen in Fig. [H As reported in Table [21 this 
AGN falls on the module#.channel# combinations 22.75, 
20.71, 4.11, and 6.15. Of these, 22.75 and 4.11 are known to 
be moderately affected by the rolling-band and moire arti¬ 
facts, while 20.71 and 6.15 are thought to be free from these 
artifacts {Kepler Instrumentation Handbook). Yet the same 
oscillatory signal is observed even during the quarters when 
the target falls on the un-affected module#.channel# com¬ 


binations, strongly suggesting that the variations are intrin¬ 
sic to the AGN. Further confirmation comes from the fact 
that kplr010841941 shares the same module#.channel# se¬ 
quence as kplr002694186 but lags behind kplr002694186 by 
one quarter, suggesting that moire and rolling-band features 
seen in kplr002694186 should also be seen in kplr010841941 
but lagging behind kplr002694186 by one quarter. No such 
phenomenon is observed, suggesting that the ripples are 
intrinsic to the AGN light curve rather than caused by 
instrumentation. In the case of kplr003337670, the mod¬ 
ule#.channel# sequence is 16.53, 24.81, 10.29, and 2.1. Of 
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Figure 4. PDCSAP light curve of the Seyfert 1 kplr002694186 
— Top panel: unstitched; bottom panel: stitched. The oscillatory 
features in this light curve may at first be suspected to originate 
in the known Kepler ‘moire’ and ‘rolling-band’ instrumentation 
effects. This AGN lands on the repeating module#:.channel# pat¬ 
tern 22.75, 20.71, 4.11, and 6.15. While 22.75 and 4.11 are known 
to be prone to the moire and rolling-band problems, the other 
two are thought to be free of these effects suggesting that these 
oscillatory features are intrinsic to the light curve itself. 


these only 24.81 suffers from moderate amounts of rolling- 
band and moire while the other 3 channels are clean. We 
observe (Fig. [2| that the rapid low amplitude oscillations 
are present in the target during the first observed quarter 
when the target fell on the module#.channel# combina¬ 
tion 16.53, which is thought to be clean. During the next 
quarter, on 24.81-a dirty module#.channel#-the AGN in¬ 
stead exhibits a lower-frequency, larger amplitude type oscil¬ 
lation. However, in the third observed quarter kplr003337670 
falls on 10.29, a clean module#.channel#, while contin¬ 
uing to exhibit the same type of behavior. Such behav¬ 
ior cannot be ascribed to rolling-band or moire effects. 
Furthermore kplr003337670 shares module#.channel# se¬ 
quence with kplr01215890 and kplr012208602, both of which 
show practically no rippling features in their light curve. The 
triplet of kplr007610713, kplr009650715, & kplr006595745 
all share the same sequence of module#.channel# but have 
very different looking light curves, suggesting that the moire 
and rolling-band effects are weak at best. Although these 
arguments do not conclusively prove that moire and rolling- 
band effects are not partly responsible for the ripple features 
observed in these light curves, we shall assume for the time 
being that these oscillatory features are real and postpone 
a more thorough investigation of residual data artifacts to 
the re-calibration of the light curves suggested in SecO 

The short-period low amplitude oscillations seen in cat¬ 
egory 3 AGN are punctuated by periods when the AGN 
appears to display a rather different sort of behavior that 



Figure 5. ‘Stitched’ PDCSAP light curve of the FSRQ 
kplr011606854. After a featureless initial period, a burst of semi- 
oscillatory behavior occurs at the ~ 75 d mark culminating in the 
flare at the ~ 200 d mark. The light curve exhibits more oscilla¬ 
tory behavior after the flare but ultimately appears to settle back 
into an almost featureless phase by the ~ 300 d mark suggesting 
that at least some AGN may display ‘on’ and ‘off’ phases in their 
light curve. 


also has an oscillatory nature but is usually of larger am¬ 
plitude and longer period, as described in the case of 
kplr003337670. Occasionally these punctuating periods re¬ 
sult in overall changes in the flux level of the object. We 
caution that these oscillatory features may not be genuinely 
oscillatory, i.e. there are classes of ARMA random pro¬ 
cesses that can generate superficially oscillatory behavior 
gi ven appropriate ARMA pa r amete r choices-see examples 
in [Woodward. Gray, fc Elliotl ll2012ll . 

Gategory 4 AGN (kplr006690887, kplr006595745, 

kplr005597763, kplr0011606854, and kplr010663134) exhibit 
flares in their fight curves as seen in row 4 of Fig. [31 Three 
out of the five members of this class are blazars, implying 
that this type of behavior may be characteristic of blazars. 
However, no flares can be observed in the 4‘^ blazar in our 
sample, kplr011021406 (which is also one of the objects that 
exhibits no variability). A solution to this puzzle may lie 
in the fight curve of kplr011606854 in Fig. [S] We see that 
after a featureless initial period, there is a burst of semi- 
oscillatory behavior at the ~ 75 d mark culminating in the 
flare at the ~ 200 d mark. The light curve exhibits more 
oscillatory behavior after the flare, but ultimately appears 
to settle back into an almost featureless phase by the ~ 300 
d mark, suggesting that at least some AGN may display ‘on’ 
and ‘off’ phases in their light curve. 

As mentioned earlier, not all of the objects show a mea¬ 
surable variability signal. Fig [6] shows stitched PDCSAP 
light curve of the Seyfert 2 galaxy kplr008703536. The vari¬ 
ations seen in this light curve appear to be purely noise. 
Similar behavior is seen in the stitched light curve of the 
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Figure 6. ‘Stitched’ PDCSAP light curve of the Sy2 galaxy 
kplr008703536. We detect no variability signal in this object - 
the observed variations in this light curve may be ascribed to in¬ 
strument noise suggesting that no endemic calibration problems 
exist in the dataset. 


FSRQ/Sy 1.5 kplr011021406. The absence of features in 
these PDCSAP light curves and in a large number of stellar 
light curves available on the MAST suggests that the latest 
PDC module (Data Release 21 onwards) corrects most of 
the systematic errors in the Kepler SAP light curves, leav¬ 
ing behind no endemic spacecraft event related artifacts. 
The only remaining sources of non-physical variability such 
as moire patterns, etc. must therefore be restricted to spe¬ 
cific pixels and CCDs, allowing us to assume with high con¬ 
fidence that our stitched light curves are primarily astro- 
physical signal. The lack of variability observed in these two 
light curves confir ms that not all AGN exhibit variability 
dSesar et al.ll200^ . though our sample is too small to be 
able to put constraints on the fraction of AGN that exhibit 
variability. The QSO kplr012208602 exhibits a very weak 
variability signal, occasionally displaying weak, intermittent 
semi-sinusoidal changes in flux. We posit that these changes 
would be un-observable using ground-based studies due to 
the weakness of the variability signal. The light curves for 
all 3 objects are plotted in the last row of Fig. |3] 

We caution that per-channel/per-pixel effects may be 
responsible for some of the features observed in our light 
curves, particularly amongst the category 2 light curves. 
A true determination of the reality of some of these fea¬ 
tures will have to await the very systematic, customized 
per-object calibration discussed at the beginning of Sec. [5] 
In the next section we discuss our method for analyzing the 
multi-quarter light curves. 


7 STRUCTURE FUNCTIONS 

We probe the variability properties of the Kepler AGN light 
curves for consistency with the AR(1) process of Sec. [2] by 


using structure functions to determine best-fit values of the 
parameter 7 in the damped power-law (DPL) model of Eq. 
®. The DPL model is used to generate mock light curves 
with the same cadence and missing observation properties 
as the AGN under investigation. We estimate best fit model 
parameters by comparing the structure function of the ob¬ 
served light curve to the ensemble of structure functions 
computed from the mock light curves consistent with the 
model parameters. 

Many definitions of structure functi ons have ap¬ 
peared in astronomy literature (see iGraham et al.l 
(|2014| for an overview). We use the d ef initio n pro- 
vided by ISimonetti, Cordes. fc HeeschenI (Il985l l and 
iRvtov. Kravtsov, fc Tatarskiil i 198'^ The n-th Order 
Structure Function of the process F{t) is 


SF^^\At) = {E^J^>{t,Atr)t, 




(15) 


where At) is the N-th increment of the process F{t) 

at time-lag At. The N-th increment is given by 

H^"^)(t,At)= (-!)"*( ^ )F(t+[iV-m]At). (16) 


Using Eq. (fTSll and Eq. (IT 6 ll . the 1®* increment and 1®* Order 
Structure Function of the process F{t) (hereafter referred to 
as ‘the increments’ and ‘the SF’) may be estimated using 


H(t,At) =F(t + At)-F(t), (17) 

and 

SF{At) = ([F{t + At)-F{t)f). (18) 


Structure functions offer benefits over directly estimating 
the PSD and ACVF/AGF of the process. Unlike the PSD, 
structure functions may be estimated in real space as op¬ 
posed to Fourier space, making them more robust estimators 
of model parameters than PSD e stimators that suffer from 
windowing and aliasing concerns llRutmanlll978li . Structure 
functions offer an advantage over estimating the ACVF and 
AGF because of the de-trending properties of structure func¬ 
tions: an n-th Order Structure Function is insensitive to (n- 
l)-th order trends in the dataset. Unfortunately, the use¬ 
fulness of this property is limited to the lower order struc¬ 
ture functions. The Kepler light curves are not long enough 
for computation of significantly higher order structure func¬ 
tions. Structure functions of all orders are related to the 
ACVF by simple linear equations. For example, the first or¬ 
der structure function can be related to the ACVF of F{t) 
by 


SF{At) = 2ACVF{0) - 2ACVF{At). (19) 


To understand what information the structure function 
conveys, consider the histograms of flux differences, i.e. T** 
increments for various values of time-lag At shown in Fig. [3 
By definition, the structure function at time-lag At is the av¬ 
erage of the squares of the values entering each histogram, 
i.e. it is the variance of the flux differences at the chosen 
time-lag At. The structure function therefore quantifies how 
the variance of the flux differences changes as a function of 
time-lag. In the case of objects that show stochastic varia¬ 
tions in F{t), the histograms look like bell-shaped curves at 
small time-lags. As the time-lag increases, the width of the 
bell-curves increase until a critical time-lag r is reached. In 
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Figure 7. Histogram plots of the increment {F{t-\-At) — F{t)) 
for various choices of the time-lag At drawn from the light curve 
of the Seyfert 1 galaxy kplr006932990. Each histogram shows the 
distribution of the increment values for a given choice of time 
lag At. The I'^^-order structure function is the variance of the 
distribution in each histogram and quantifies how the variance of 
the increments changes as a function of time lag. 


the case of the AR(1) process, this critical time-lag can be 
identified with the the turnover timescale in Eq. (Ell- 

Astronomical time series have measurement noise 
that is usually modeled as uncorrelated white-noise, i.e. 
ACVF{0) = P + (j% where cr^ is the contribution to the 
variance of the time series at zero-lag from the actual signal, 
while characterizes the noise level of the measurement 
noise. In the case of the AR(1) process, the ACVF takes the 
form in Eq. 0. If the variability observed in Kepler AGN 
light curves is well described by an AR(1) process, then the 
structure function should be well described by 

SFDRw{At) = 2P{l-e-\^\) + 2a%. (20) 

This may be tested by estimating the structure function of a 
real AGN light curve to check if it is consistent with Eq. (1201) 
or with a more general form based on the DPL PSD of Eq. 
®. Unfortunately there is no closed-form expression for the 
auto-correlation function corresponding to the DPL PSD, 
making it impossible to directly fit the functional form for 
the structure function observed for a given AGN. In the next 
section, we describe a Monte-Garlo estimation technique for 
recovering the DPL model parameters. 


8 MONTE-CARLO ESTIMATION OF THE 
DAMPED POWER-LAW MODEL 
PARAMETERS 

We estimate the DPL parameters in Eq. via Monte-Carlo 
simulations. We compare the real structure function of an 
AGN to simulations computed from mock light curves gen¬ 
erated using the DPL model of Eg. ll^. We gen erate ‘mock’ 
light curves using the iTimmer fc (Il995l l method. To 

create a single mock light curve, pseudo-random numbers 
are generated using the Fast Mersenne Twister SFMT19937 
generator seeded with hardware-generated random numbers 
(generated using Intel RDRAND instruction) to ensure that 



Figure 8 . Observed light curve of the Sy 1 AGN kplr006932990 
(orange) along with an example mock light curve (light green) 
generated as described in[ 8 ] (generated with 7 = 2.7 and r = 27.5 
d). The mock light curves used in Monte-Carlo estimation pro¬ 
cess have the exact sampling pattern of the observed light curve 
ensuring that the structure functions of the mocks are subject to 
the same sampling issues as that of the real galaxy. 


the ran dom number seque nces are free of artificial corre¬ 
lations (iGoddingtonI Il99-d) induced by poor random seed 
choices. At this intermediate stage, the mock light curve is 
oversampled by a factor of 10 i.e. we generate points at 10 x 
the required cadence in order to avoid sampling issues. To 
include low frequency modes that are not adequately char¬ 
acterized by the length of the observed light curve, the inter¬ 
mediate mock light curve is much longer than is ultimately 
required to make the final mock; mock light curves gener¬ 
ated in this manner are capable of exhibiting low frequency 
modes longer than the length of the observed data. FFTs 
are most efficient for data sequences that are a power of 2; 
for this reason, we pick the intermediate (including the over- 
sampling) to be of length 2^®. This results in the intermedi¬ 
ate mock light curve being between 15 x to 45 x the length 
required for the final mock light curve depending on the ac¬ 
tual length of the observed light curve. We pick a uniformly 
distributed random segment of the intermediate overly-long 
light curve that has the same length as the observed light 
curve and generate another stream of un-correlated Gaus¬ 
sian random deviates to simulate the white-noise properties 
of Kepler instrumentation noise. After adding this ‘measure¬ 
ment noise’, we set data points corresponding to the un¬ 
observed cadences in the real light curve to 0. This proce¬ 
dure creates a final mock light curve with identical sampling 
and noise properties to the real light curve. Fig.[8]shows the 
true light curve (orange) along with an example mock light 
curve (light green) for the Sy 1 AGN kplr006932990 illus¬ 
trating what the mock light curves look like for the best fit 
DPL parameters for this object. 

We compute the structure functions using Eq. m 
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Figure 9. Distribution of mock SF{St) estimates for 1000 simu¬ 
lated light curves of length T = 10 d, with 7 = 2.0 S>c r = 7.5 d for 
= 5 d. The upper panel shows the (clearly non-Gaussian) dis¬ 
tribution of SF{5t), while the lower panel shows the distribution 
of \ogiQ SF{St). The Shapiro-Wilk test W- Sz p-values confirms 
that the distribution of log;i^Q SF{5t) is closer to Gaussian than 
that of the raw SF{St). 


Figure 10. Observed structure function of the Sy 1 AGN 
kplr006932990 (orange), best-fit model structure function with 
standard deviations (purple) and the mock structure function 
(light green) corresponding to the mock light curve in Fig. [S] 
The value of 7 is most sensitive to the short timescale (At < 5 d) 
behavior of the structure function where the variance allowed is 
small. 


modified to account for missing values - 


SFn 


E N — n 
i=l 


WiWi+n 


[Fi + n — 


E N — n 
i=l 


WiWi+„ 


( 21 ) 


where n = At/St is the lag in cadences as opposed to physi¬ 
cal time. The weights are defined to be iCi = 1 for observed 
cadences and Wi = 0 otherwise. 

It is known that, for any given lag At, the distribn- 
tion of SF{At) gen erated nsing the Timmer-Konig method 
is n ot Gaussian dEmmanoulopoulos, McHardv, fc UttlevI 
I2OIOI I. However, the logarithms of the structure function, 
log,Q SF{At), are distributed as per a Gaussian distribution. 
As can be seen in Fig. [9] histograms of mock SF{At) and 
log,Q SF{At) values for a set of 10000 mock light curves con¬ 
structed with PSD slope 7 = 2, and characteristic timescale 
T = 7.5d at At = 5d for mock light curves of length T = 10 d 
suggest that while the estimates of SF{At) are not Gaussian 
distributed, the estimates of logjQ SF{At) are. This obser¬ 
vation may be confirmed using a test such as the Shapiro- 
Wilk test. Note however, that even though the Shapiro-Wilk 
test confirms that logjQ SF{At) is closer to Gaussian than 
SF{At), it is not always Gaussian in an absolute sense; a 
more through investigation of the distribution properties of 
the structure function may be warranted. We use logjQ SFn 
instead of SFn to compare structure functions between ob¬ 
servations and simulations. 

Since the logarithms of the structure function estimates 
are drawn from Gaussian distributions, we can use the mean 
and standard deviation vectors of the logarithms of the mock 
structure functions to estimate the best-fit DPL model pa¬ 
rameters. We compute the mean Hn = (logjQ SFn) and vari¬ 


ance cr^ = {{logiQ SFn — {logiQ SFn))'^) of the ensemble of 
the logarithms of the mock structure functions and numeri¬ 
cally minimize 

^2 _ inn ~ logic SFn^obs) ^^2'^ 

n=l 

using the Gonstrained Optimization BY Li near Approxi - 
mations (COBYLA) optimization algorithm (IPowelll [l994li 
to search the parameter space for the best-fit estimates of 
7, r, as, and crjv. Coarse optimization is performed using 
10 "^ mocks at each step in the optimization, after which the 
best-fit parameter estimates are refined using 10® mocks at 
each step in the optimization process. Such large numbers 
of mocks are required to ensure that values computed at 
the same point in parameter space (using different choices 
of random seeds) are within the tolerances chosen in the 
optimization step. 

The resulting values of the reduced chi-square per 
degree-of-freedom, Xdof are close to 1. We also compute 
the percentile, P, of the value of XooF corresponding to the 
best-fit parameter values using a fresh set of 1000 mock light 
curves. Percentile values P close to 100% mean that the val¬ 
ues of Xmock are almost uniformly smaller than Xobs for the 
best-fit model parameters, indicating that the model is a 
poor fit to the data. 

Fig. [TO] shows the observed structure function of the 
Seyfert 1 kplr006932990 (orange) along with the best-fit 
mock structure function and standard deviation (purple) 
and the mock structure function (light green) correspond¬ 
ing to the mock light curve in Fig. [8| On short timescales 
(At < 5 d), the mock structure functions show very small 
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standard deviation, while on longer timescales, the structure 
function is allowed to vary more significantly. The value of 
7 is determined almost exclusively by this short timescale 
(At < 5 d) behavior of the structure function where the 
variations allowed in the possible structure function realiza¬ 
tions are small. Note that in the case of light curves that 
have 7 > 2.0, as the PSD turns over from shorter timescales 
(At ^ 5 d with local PSD slope 7 > 2) to longer timescales 
(At ^ 10 d with local PSD slope 7 —>■ 0), the local value 
of the PSD slope passes through 7 = 2. Most ground-based 
surveys, such as the SDSS, have large (compared to Kepler) 
photometric errors and are measuring the structure function 
on timescales At ^ 5 d resulting in estimates of 7 2. 

Ideally one would estimate confidence intervals for the 
model parameters using algorithms such as MCMC to sam¬ 
ple the parameter. Unfortunately, the Monte Carlo simula¬ 
tion is computationally very expensive with the calculation 
of taking 3 h for every location in parameter space on 
a 16-core Xeon machine with 2 Xeon Phi accelerator cards. 
This makes the Monte-Carlo technique of fitting the struc¬ 
ture function unsuitable for an MCMC style analysis. To 
estimate the error in our computation of the PSD slope 7 
we generated simulations of short segments of mock light 
curves with known DPL model parameters and then pro¬ 
ceeded to recover the input parameters. Based on this, we 
estimate that the error in the value of 7 obtained is on the 
order of 10%. 


9 DETERMINATION OF Tmin 

We estimate the shortest timescales on which we observe 
variability in the light curves. Our estimation of the DPL 
model parameters in Sec. [8] yields an estimate of the un¬ 
correlated white-noise level present in each light curve. This 
noise level can be seen in the structure function in the form 
of a flattening-out of the structure function on very short 
timescales. For example, in Fig. 1111 we observe that the 
structure function of the Seyfert 1 kplr006932990 begins to 
level off to a value of 100[e~]^ on timescales of 0.5 d. 
From Eq. (PJ)). we expect that this is caused by the noise 
level i.e. 2a% 100[e“]^. By removing this noise level from 

the observed structure function, we can estimate a noise¬ 
less version of the structure function. To estimate the short¬ 
est timescale on which we observe variability in the light 
curve, we find the time-lag at which this noiseless structure 
function crosses the noise floor. In Fig. [TT] we plot the ob¬ 
served noisy structure function of kplr006932990 (orange), 
our best fit estimates of this structure function (purple), 
and the noise level determined from the structure function 
fit (red dashed-line). After removing this noise level from the 
observed structure function, we obtain the noiseless struc¬ 
ture function (red). We also show the mock structure func¬ 
tion (light green) corresponding to the mock light curve in 
Fig. 0 along with the noiseless version (green) of this mock 
structure function. We define the variability onset timescale 
Tmin to be the intersection point (indicated by blue dashed- 
line) of the noiseless structure function with the noise floor. 
We compute the value of the variability onset timescale for 
the AGN in this study along with the rest-frame sampling 
interval for each object and present the results in Table (3] 



logio^hest (d) 

Figure 11. Observed structure function of the Sy 1 AGN 
kplr006932990 (orange) and the best-fit model structure function 
with standard deviation (purple). The noiseless structure function 
(red) is obtained by removing the noise floor (red dashed-line) 
after which the shortest variability timescale Tmin (blue) is de¬ 
termined from the intersection of the noiseless structure function 
with the noise floor. Tmin typically < 1 d implies that future AGN 
variability studies should attempt to sample the short timescale 
variability properties of AGN. 


10 OBSERVED STRUCTURE FUNCTIONS 
AND FITS 

The excellent short timescale sampling properties of Kepler 
make it possible to study the AGN structure functions with 
unprecedented levels of detail. Using the Monte-Carlo fitting 
method of Sec. [S] we fit the DPL parameters of Eq. (0 for 
the 17 most variable objects. Fig. [12] shows the observed 
structure functions (orange), Monte-Carlo hts (purple line), 
and la variation (purple shaded region). 

Sampling information and estimates of the DPL param¬ 
eters are given in Table O for each Kepler AGN. Column 1 
lists the KeplerlD of the object. Column 2 lists the length 
of each light curve in quarters. Each quarter spans about 3 
months in duration and has ~ 4300 data points. Included 
in this hgure are quarters during which the target fell on a 
failed CCD module, resulting in one or more missing quar¬ 
ters of data if flux measurements exist before and after the 
‘missing’ quarter. Column 3 lists the rest-frame sampling in¬ 
terval for each object in minutes. The next two columns list 
the best-fit values for the DPL parameters 7 and r from Eq. 

where r is measured in rest-frame days. Column 6 lists 
the best-fit estimate for the minimum timescale on which 
variability is observed Tmin in minutes. The ‘goodness-of- 
ht’, i.e. the reduced chi-square, is listed in Column 7 while 
Column 8 lists the percentile of the chi-square of the ob¬ 
served structure function fit as compared to simulations 
drawn from the DPL model i.e. it is also an indicator of 
‘goodness-of-ht’ in that a very large value of the chi-square 
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Figure 12. Observed structure functions (orange) and fits using mock light curves (purple) of all the AGN. Note that all the structure 
functions have been plotted at the same scale. 


percentile indicates that the DPL model is insufficient to 
explain the data. 

We observe that in all but one case (the blazar 
kplr00669887), Kepler samples the fight curve on shorter 
sampling intervals direst (Col. 4 in Tab. [3]) than Tmin (Col. 
7). The median variability onset timescale is ~ 23 h in the 
rest frame of the object, while the minimum and maximum 
values are ~ 2.5 h and ~ 1.6 d respectively. The structure 
function of kplr006690887 is notable for being bereft of a 
well-defined noise floor implying that this BL Lac object 
varies even on the 25 min timescale probed by Kepler in the 


rest-frame of this object. This suggests that the Kepler sam¬ 
pling pattern is more than adequate to study optical AGN 
variability. We caution that Tmin is merely an upper esti¬ 
mate for the shortest timescale on which variability can be 
observed given the noise level, therefore even better photo¬ 
metric precision than is available with Kepler is required to 
find the true shortest timescale on which AGN vary. 

We categorize the fight curves of our AGN based on the 
presence of oscillatory features and flares in Sec.[Sl All 3 ob¬ 
jects that we grouped in category 1 (completely stochastic¬ 
looking fight curves-row 1 of Fig. El are fit by DPL models 
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Table 3. Sampling Parameters Sz Results. 


Identifier 

KeplerlD 

Sampl. Param. 

L direst 

(#Q) (min) 

7 

Analysis Results 

Train X^oF 
(d) (min) 

P 

% 

6932990 

14 

28.64 

2.7 

27.5 

163 

1.43 

83.3 

12158940 

12 

27.58 

2.5 

40.2 

440 

1.15 

71.3 

11178007 

11 

27.27 

2.9 

16.8 

428 

0.75 

38.5 

9650715 

4 

26.11 

2.2 

47.8 

502 

0.64 

40.2 

2837332 

7 

26.04 

2.1 

41.0 

1358 

0.86 

55.4 

3347632 

7 

25.41 

2.5 

18.0 

1453 

1.07 

68.1 

5781475 

4 

24.08 

3.1 

81.0 

1619 

1.18 

75.3 

2694186* 

11 

27.02 

2.1 

134.3 

2339 

1.72 

88.6 

9215110 

8 

26.39 

2.2 

55.1 

1704 

3.22 

98.6 

10841941 

7 

25.54 

2.3 

10.9 

2151 

1.23 

78.3 

3337670 

7 

21.51 

2.2 

125.4 

9946 

0.16 

1.80 

7610713 

8 

20.45 

2.5 

29.3 

1677 

0.63 

35.2 

6690887 

7 

25.50 

2.2 

3.1 

- 

13.37 

99.9 

6595745 

7 

19.59 

1.7 

76.8 

711 

1.23 

76.9 

5597763 

7 

18.11 

1.8 

12.5 

1151 

3.14 

98.7 

11606854 

12 

15.28 

3.0 

1.6 

2005 

8.85 

99.9 

10663134 

12 

11.68 

2.1 

21.7 

463 

5.47 

99.9 

8703536 

12 

27.92 

- 

- 

- 

- 

- 

11021406 

12 

19.45 

- 

- 

- 

- 

- 

12208602 

12 

14.02 

- 

- 

- 

- 

- 


*kplr002694186 has a total light curve spanning 17 quarters in 
the MAST database. Data is available only for quarter 0 (as an 
exo-planet search program target) and then subsequently from 
quarters 7 through 17 (as a GO program target). We discard the 
quarter 0 light curve segment because of the unreliability of the 
photometry during this initial phase of the mission and report 
kplr2694186 as having a light curve of length 11 quarters. 

The AGN are grouped based on light curve features as in Table 
^ Column 2 (‘ 7 ^ Q’) is the approximate length of each light curve 
in quarters. Note that occasionally a target will land on one of 
the failed CCD modules and will not have data available for the 
quarters corresponding to such occurrences. Column 2 includes 
such ‘missing’ quarters since the overall length of the light curve 
just depends on the first and last quarter observed. For exam¬ 
ple, the radio-loud AGN kplr011021406 does not have data for 
quarters 8 , 12 , and 16 but since this object was observed starting 
in quarter 6 and ending in quarter 17, the light curve still has 
a total length of 12 quarters of data. Other than the incomplete 
17th quarter, each quarter contains approximately 4300 individ¬ 
ual data points making the longest light curve (kplr006932990) 
over 60,000 points long. Column 3 lists the rest-frame sampling 
interval for each object. Columns 4 through 6 list the results of 
the structure function analysis where 7 is the best fit slope of the 
power-law portion of the power spectral density, r is the turnover 
timescale beyond which individual data points are essentially un¬ 
correlated, and Train IS the shortest timescale on which variability 
is observed. Columns 7 8 list measures of the ‘goodness-of-fit’ 

i.e. the reduced ^^e percentile value P of the chi-square 

compared to chi-squares computed from 1000 mock light curves. 

with median 7 = 2.7 with xi>oF close to 1. This value of 
7 is significantly different from that expected of an AR(1) 
process and indicates that the light curves of these 3 objects 
are smoother than a damped random walk {'Xdrw = 2 . 0 ). 
Prominent in the observed structure function fits for these 
three AGN in the 1^* row of Fig. [12] is the presence of a 
‘dip’ feature on time-lags At between ~ 10 and 100 d. This 
dip feature may be interpreted as an excess of correlation on 
these timescales. The dip starts at ~ 10 d in the rest-frame 
of the objects, suggesting that it is probably intrinsic to the 
light curves of these objects. Such dips rarely occur in the 



loglO^t^est (d) 

Figure 13. Observed structure function of the Sy 1.9 AGN 
kplr009215110 (orange), noise-removed observed structure func¬ 
tion (red), example mock structure function (light green), noise¬ 
less mock structure function (green), and best-fit model structure 
function with standard deviations (purple). This structure func¬ 
tion exhibits varying PSD slope 7 . The structure function rises 
steeply when Atrest 3 d after which is flattens somewhat. The 
DPL model is incapable of reproducing such behavior. 


structure functions of the simulated light curves (generally 
these exhibit wiggling behavior after the turnover timescale 
is reached). This suggests that, even though all three objects 
have values consistent with the DPL model (P ranges be¬ 
tween 40% and 80%), the DPL model is probably too simple 
to characterize the variability observed in these objects. 

Of the objects categorized as intermediate between cat¬ 
egories 1 and 3 (stochastic-looking light curves-l-oscillatory 
features-row 2 of Fig. 01 , the first 2 have PSD slopes (7 = 2.2 
and 7 = 2.1 respectively) close to that expected from an 
AR(1) process. Of these first two objects, kplr009650715 
has a fairly low XdoF ~ 0.64 implying that the DPL model 
may be overfitting the light curve of this object. The lat¬ 
ter two objects, kplr003347632 & kplr5781475, have much 
steeper PSD slopes (7 = 2.5 and 7 = 3.1 respectively) com¬ 
pletely inconsistent with the AR(1) process and are well fit 
by the DPL model. The median PSD slope for this category 
is 7 = 2.35 and the individual structure functions are shown 
in row 2 of Fig. [12] 

Objects that fall into category 3 (strong oscillatory fea¬ 
tures in light curves-row 3 of Fig. 0 based on a visual in¬ 
spection of the light curve have the largest number of light 
curve slopes consistent with that expected from the AR(1) 
process. Four out of 5 have 7 ~ 2.2, while the median esti¬ 
mate of 7 = 2.2 for all 5 combined. Individual estimates of 7 
range between 2.1 (kplr002694186) and 2.5 (kplr007610713). 
Even though our structure function estimation algorithm 
produces estimates of the turnover timescale r for all of these 
objects, a visual inspection of the corresponding structure 
functions in the row 3 of Fig. [T^ suggests that the turnover 
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timescale is suspect in all but the case of kplr010841941, as 
borne out by the mock structure function error contours. 
Dips similar to those observed in category 1 objects are also 
seen in the structure function of some of these objects. As in 
the case of category 1 AGN, the dip feature begins around 
10 d in the rest-frame of the objects and is usually ab¬ 
sent by ~ 100 d. kplr002694186 & kplr009215110 exhibit 
structure functions that appear to have different slopes on 
different timescales. This presence of multiple slopes may 
be responsible for the relatively poor quality of the fits of 
the structure functions of these objects (xdoF is 1-72 and 
3.22 respectively). Fig. 1131 shows the structure function of 
the Sy 1.9 AGN kplr009215110. On short timescales < 3 d 
the structure function rises steeply, after which it flattens 
somewhat. The DPL model is incapable of producing such 
behavior. Similar behavior is also observed in kplr006595745. 

Category 4 objects with flares in the light curves ex¬ 
hibit a broad spread in estimated 7 ranging from a fairly 
flat 1.7 (kplr006595745) to a very steep 3.0 (kplr011606854). 
However, this category of objects have the worst fits as a 
class, with very suspicious XdoF values such as 13.4, 8.9, & 
5.5 (kplr006690887, kplrOl 1606854, & kplr010663134). All of 
these objects have much higher values than is the norm 
for mock light curves generated using the DPL model, sug¬ 
gesting that the DPL model (and therefore the AR(1) pro¬ 
cess) are wholly unsuitable for modeling the light curves of 
objects that show flaring behavior. 

The last row of Fig. [12] shows the individual 
structure functions of kplr008703536, kplr011021406, and 
kplr012208602. A visual inspection of the corresponding 
light curves in Fig. [3] suggests that these objects do not 
exhibit significant levels of variability. This observation 
is borne out by the structure functions of these objects, 
which are mostly flat and featureless for the duration of 
observed time-lags. On long timescales 100 d) the 
structure function begins to exhibit low amplitude ‘wig¬ 
gles’. It_js_well_Jcnown from the theory of ACF estima¬ 
tion llBrockwell &: Davisll2002l l that similar wiggles occur on 
timescales comparable to the length of the time-series and 
are not significant. 


11 COMPARISON 

We observe that not all AGN exhibit DPL parameter 7-^2, 
consistent with an AR(1) process. Similar observations have 
been reported by others. In this section we compare our 
results with those from other AGN variability studies per¬ 
formed using both ground- as well as as space-based instru¬ 
ments. We begin by comparing our results with other studies 
of Kepler AGN. 

11.1 Kepler AGN Variability Studies 

We have shown that the AGN in this study exhibit PSD with 
shapes that are superficially similar to that expected from a 
simple variability model, such as the damped random walk 
or AR(1) process—a flat power spectrum on long timescales 
that turns into a 1/P power-law section after frequencies 
higher than 1/r. However while the AR(1) requires the value 
of 7 to be exactly 2, we find a wide r ange of 7 values incon¬ 
sistent with the AR(1) model. Both [^rini fc Rvld (l2012l l 


and iMushotzkv et all ll201lll obtained values of the PSD 
slope for the Seyfert I AGN Zw 229-15 or kplr006932990. 
ICarini fc Rvld (l2012|j applied various PSD models and found 
that the best fit estimate of the PSD slope is 7 ~ 2.8 for 
both their ‘knee model’ as well as their ‘broken power-law’ 
model. Their estimates of the break timescale ranged from 
r = 92 ± 27/21 d for the ‘knee’ model to r = 43 ± 13/10 d 
for their ‘broken power-law’ model, which is in good agree¬ 
ment with our estimate of 7 = 2.7 while w e find a somewhat 
shorte r characteristic timescale r = 27.5 d. IMushotzkv et al.l 
1I2OIIII computed the PSD slope for 4 individual quarters 
and found an average PSD slope of 7 = 3.115, significantly 
steeper than our estimate. It should be noted that these re¬ 
sults have been obtai ned using different ve rsions of the final 
Kepler data product. ICarini fc Rvld (I2OI2I') were able to use 
ground-base d photometry from the reverberation mapping 
campaign of iBarth et al.l 1I2OIIII along with a customized 
aperture to stitch t ogether the 4 quar t ers of light curve data 
used in the study. IMushotzkv et all (1201 il l calculate PSD 
fits for 4 individual quarters of Kepler data but use the 
un-processed SAP light curve in their estimate of the PSD. 
We used PDCSAP fluxes (Data Release 21-I-) from all 14 
available quarters of data in our analysis. These fluxes, cali¬ 
brated using the newest version of the PDG module, should 
suffer only marginally from instrumentation issues. Thus, it 
is likely that the range of PSD estimates observed is a re¬ 
sult of the subtle but important differences between the data 
products used. Given that all 3 studies have used different 
analysis techniques but arrive at similarly steep PSD slope 
estimates, it is clear that the PSD slope is robustly steeper 
than 1/ P and is irreconcila ble with the AR(1) process. Most 
recently. iKellv et al.l (l2014fl attempt to use the Kalman filter 
to apply a general ARMA model to a single quarter of data 
from kplr006932990 and find that the observed variability is 
best-fit by an ARMA process with 6 Auto-Regressive and 
4 Moving Average components—a much more complicated 
stochastic process than a simple damped random walk. 

IMushotzkv et al.l 1I2OIIII also obtained PSD slope es¬ 
timates for 3 other AGN: kplr012158940, kplr011178007, 
and kplr002694186. Their analysis was performed on the 
SAP fluxes observed during individual quarters. The av¬ 
erages of their estimates of the PSD slopes: 7 = 2.67 
(kplr012158940) and 7 = 2.92 (kplr011178007) agree fairly 
well with our estimates in Table |3] however, they observe 
a much steeper 7 = 2.845 for kplr002694186 than our es¬ 
timate of 7 = 2.1. Both kplr012158940 and kplr011178007 
exhibit strong variability-most of the ‘signal’ observed in 
the uncalibrated SAP light curves for both objects was in¬ 
trinsic to the source and retained in the calibrated PCD- 
SAP light curves. On the other hand, the calibrated PDG- 
SAP light curve of kplr002694186 lacks the large amplitude 
flux variations seen in the original SAP light curve, indicat¬ 
ing that those variations were instrumental effects that may 
have inadvertently intr oduced power at l o wer fr equencies in 
the PSD computed by IMushotzkv et al.l ll201lll . artificially 
making the PSD slope steeper than that computed in our 
analysis. 

lEdelson et al.l ll2013h have examined the variability 
properties of the light curve of the blaz ar kplr0066908 8 7. Us - 
ing segments of data spanning ~ 5.5 d. lEdelson et al.l (l2013fl 
found that while power law PSDs yield unacceptable fits for 
the observed PSD, the least unacceptable model is a bending 
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power law PSD with slopes 7 = 1.30±0.14 and 2.86±0.34 if 
the highest frequencies are excluded from the analysis. Us¬ 
ing the full light curve, we find that our fit is of very poor 
quality with 7 = 2.2. However, flaring is not an AR(1) pro¬ 
cess. An AR(1) process is completely incapable of produc¬ 
ing the flare features observed in t his AGN. Ground based 
studies of blazar variability, such as iRuan et al.l (|2012| ~I. have 
attempted to model blazar variability with an AR(1) pro¬ 
cess using 101 blazar light curves drawn from the LINEAR 
near-Earth asteroid survey. The LINEAR survey collected 
time series data over a period of 5.5 yr with two 1.01 meter 
telescopes giving a r-band survey depth of 18 mag at 5 a. 
Sources in the survey within 10° of ecliptic have 460 obser¬ 
vations on average, while sources further from the ecliptic 
only have 200 observations over the duration of the survey, 
although the caden ce can occas i onally be as high as once 
every 15 min. While [Ruan et ahl (120121 ') find that the blazar 
light curves observed by them are well fit by the AR(1) 
process with decorrelation timescales ranging from < 1 d 
to ~ 1000 d with a peak at 100 d, we caution that the 
LINEAR survey suffers from both highly irregular, sparse 
sampling patterns, and large photometric uncertainties. 

IWehrle et al.l (l2013l 'l and iRevalski et al.l (l2014l l 
have estimated the PSD slopes of the radio-loud 


AGN kplr011021406, kplr 011606854, kplr012 28602, and 
kplr010663134 in Table [S] IWehrle et al.l ll2013l ) perform a 
PSD analysis of the un-corrected SAP flux light curve of 
each object on a per quarter basis. They also supply the 
same analysis for ‘corrected’ and ‘end-matched’ versions 
of the light curve for each quarter, where they corrected 
the SAP data by removing some of the instrumental effects 
and windowed the data by removing a linear trend from 
the data to make the first and last point of each quarter 
have identical flux value. In the case of kplr011021406, 
they concluded that this Seyfert 1.5/FSRQ is variable and 
derive a mean PSD slope 7 = 1.8 with standard deviation 
calculated from each quarter equalling 0.2. Individual 
estimates of their PSD slopes for each quarter range from 
1.8 on the low side to 2.0 on the higher end. In comparison, 
using the PDCSAP light curve from Data Release 21-|-, we 
find no intrinsic variability in the light curve of this object 
suggesting that it is essential to use the PDCSAP flux when 
studying light curve variations in order to properly remove 
spacecraft induced trends from the data. For the FSRQ 
kplrOl 1606854, they derive a mean PSD slope of 7 = 2.0, 
which is significantly flatter than the value estimated by 
us (7 = 3.0). The remaining AGN, kplr010663134, has 
slope 7 = 1.9 which in g ood agreement with o ur result 
in Table (3] Subsequently IRevalski et al.l (l2014l ~l used a 
stitching algorithm in their analysis of the multi-quarter 
light curves of these AGN. After stitching toge ther 11 
quarters of SAP flux data, IRevalski et al.l ll2014l ') obtain 
PSD slope estimates of 7 = 1.9 (kplr011606854) and 7 = 1.7 
(kplr010663134)-significantly different from the values in 
Table [S] This discrep ancy is partial l y cau sed by the use 
of the SAP fluxes by IRevalski et al.l (l2014h as opposed to 
PDCSAP fluxes by us. More importantly, the DPL model is 
just unable to fit the observed structure functions of AGN 
that have pronounced flares in the light curve. 

In conclusion, we find that our results are in good agree ¬ 
ment with the previo u s studies of Mushotzk_y_gt__aI, QOllf); 
lEdelson et al.l J2013tl : IWehrle et al.l ( 2013l ~): IRevalski et ^ 


(l2014h with most differences being attributable to the vari¬ 
ations in data-length and data-versions between the light 
curves used by us and previous efforts. 


11.2 Ground-Based Studies that Use the AR(1) 
Process 

Previous variability studies using ground-based data suggest 
that the AR(1) process is adequate to characterize AGN 
variability given the measurement uncertainties and sam¬ 
pling pattern typical of ground-based studies. The original 
estima tion performed bv iKellv. Bechtold. fc Siemiginowskal 
(I2OO9II using th e standard state-space re presentation of the 
AR(1) process dBrockwell &: Davijl2002l indicated that the 
R- and B-band optical variability observed in their sam¬ 
ple of 109 quasars and Seyferts was well modeled by AR(1) 
processes with a range of decorrelation timescales strongly 
indicative of a connection between thermal fluctuations oc¬ 
curring on the AGN accretion disk and the observed vari¬ 
ability. The 109 AGN used in this study include 59 MAGHO 
quasars observed with sampling rates of once every 2-10 d 
over 7.5 yr, 42 PG quasars observed once every 40 d over 
7 yr, and 8 Seyferts observed once every 1-10 d resulting 
in their data set having no cadence higher than ~ 1 d“^. 
The best-fit characteristic timescales for their objects range 
from -^5-20000 d with a median value of 540 d. It should be 
noted that while the shortest interval between 2 flux mea¬ 
surements in the dataset may be ~1 d, as is the case in all 
ground-based studies, this interval is not constant, i.e. the 
median sampling interval is substantially higher than ~ 1 
d. Looking at the Kepler structure functions in Fig. 1121 this 
dataset and other ground-based datasets do not sample the 
structure function very well below ~ 10 d making these stud¬ 
ies less sensitive to the short timescale (< 5 d) properties of 
AG N light curves. _ 

iKozlowski et al.l (l2010l) used the Press-Rybicki-Hew itt 
maximum likelihood technique (iRvbicki fc Pres^ Il992h to 
model a much larger sample of ~ 88700 I- and V-band candi¬ 
date quasar light curves from the OGLE-II and -III surveys 
as a set of multivariate Gaussian distributions parametrized 
by correlation matrices with the structure expected from 
the AR(1) process. The OGLE surveys that the data were 
obtained from collectively span about 12 yr and offer be¬ 
tween 100 (V-band) to 750 epochs per light curve, giving an 
average sampling rate of once every 6 d in the I-ban d. By 
using a set of known quasars IKozlowski et al.l (I2OI0I I were 
able to suggest a parameter-space cut to select quasars. The 
characteristic timescales found by them range from ~10- 
3000 d and are in agreement with the range observed by 
iKellv. Bechtold. fc Siemiginowsl^ (l2009l L 

Similar results were obtained using data from the SDSS 
Stripe 82 data set. The SDSS Stripe 82 data set spans about 
10 yr and provides about 65 epochs of data for each ob¬ 
ject in the survey giving a sampling rate of about 6-10 
observations every year. When comparing results from the 
SDSS survey to our work, it should be noted that the typi¬ 
cal SDSS quasar is much more luminous (Mj < —24) than 
the AGN observed by Kepler. iMacLeod et al.l ll2010li used 
the PRH algorithm to model the light curves of the 9000 
quasars in the SDSS Stripe 82 time domain survey as AR(1) 
processes an d observed the same range of characte ristic 
timescales as iKellv. Bechtold. fc Siemiginowsl^ 1I2OO9II and 
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iKozlowski et alJ (l 2010 l ): most quasars have decorrelation 
timescales in the range of 3 — 3000 d. These results were 
put to the test by combining data from the POSS survey 
50 yr baseline with a handful of points per light curve) 
with the Stripe 82 data bv iMacLeod et ^ ll2012ll . These au¬ 
thors found that the AR(1) process continued to fit the ob¬ 
served light curves well and yielded decorrelation timescales 
in the r ange of 5 — 2000 d for t h e light curves. Ba s ed on these 
results. iM acLeo d et alJ ll201lll , iButler fc BloomI (1201 il l , and 
IChoi et~ ( 2014) have proposed that quasars may be se¬ 
lected from variable point sources, such as non-periodic vari¬ 
able sta rs, using the vari ability parameters as a selection 
criteria. iKellv et al.l (l20l4 have shown how variability may 
be used as an estimator of the black-hole mass of the AGN. 
The AR(1) process has even been successfully applied by 
ISobolewska et aP (l20l4 to 7 -ray AGN variability observed 

by the Fermi space telescope . _ 

_ Recently Z^jt_^ljj20^) , Andrae. Kim, fc Bailer-Jonej 

(|2014 . and iGraham et al.l (I20l4 have rigorously tested 
the AR(1) process on the timescales probed by the sur¬ 
veys used in the previous studies . Using the PRH formal¬ 
ism adopted by IKozlowski et al.l (l 201 of l. but using more 
generalized for ms for th e struc ture of the covariance ma¬ 
trix employed, IZu et all (I20l4 concluded that while the 
OGLE light curves used in the study are consistent with 
the AR(1) process on the longest timescales probed by the 
OGLE data 7 yr), on short timescales (less than a few 
months) there is some evidence for a steeper than 1 //^ 
PSD. Furthermore, they concluded that while the AR(1) 
process produces acceptable fits to the light curves, on the 
whole the simple models used in that study ar e not well con¬ 
strained by the light curves. IZu et all (I20l4 concede that 
the quality of the photometry may be responsible for the 
mixed results obtained by them when fitting the various co- 
variance matrix models, i.e the scatter observed by them in 
the slope parameter fit may be either intrinsic to the light 
curves and evidence for varied behavior between the ob¬ 
jects, or a consequence of errors in the determination of the 
photometric uncertainties. A comprehensive and systematic 
Bayesian s tudy of SPSS Stripe 82 quasa r light curves per¬ 
formed by lAndrae. Kim, fc Bailer-JonesI (l2013l 'l concluded 
that more sophisticated models, such as AR(2), GARCH(l), 
and ARMA(1,1), are not required to model light curves. Of 
the 6304 quasars examined in this study, an overwhelming 
majority-5047 light curves-were best modeled by an AR(1) 
process while only a handful required a more sophisticated 
approach. 

While a large number of studies have indicated that 
the AR(1) process is sufficient to model AGN variability, 
these studies have used comparatively sparsely sampled light 
curves (10 — 100 points over several years) that do not well 
sample the light curve on short timescales (At < 5 d). 
There are indications that more complex stochastic models 
are required to chara cterize variability on short timescales. 
ISchmidt et al.l (l 2012 l ~l use the simple power-law form of the 
structure function given in Eq. (HD) to model the structure 
function of Stripe 82 quasars. They compute the structure 
function using 


SF{AU,j) = -bcr2)At, (23) 

where Am is the magnitude difference over the time interval 


At. While our structure funct ions is quadratic in fl ux differ¬ 
ence, the structure function of ISchmidt et al.l (I2OI2II is linear 
in magnitude difference. Our correlation strength parameter 
7 appears in the light curve PSD model, while th e correla¬ 
tion strength parameter A of ISchmidt et ah] (I2OI2I ') appears 
in the structure function model making it i mpossible to di¬ 
rectly compare our 7 with their A. However, I Schmidt et al.l 
(I2OI2II find that Stripe 82 quasars are fit by A values rang¬ 
ing from 0.0 to 1.2 with a peak at A = 0.25. The range 
of A values found bv ISchmidt et aP (l2012h suggests that a 
simple stochastic process model li ke the AR(1) pr o cess i s un¬ 
likely to work well for all objects. iGraham et al.l (I20l4 use 
a wavelet analysis technique known as the Slepian Wavelet 
Variance (SWV) to study 18000 quasars from the Catalina 
Real-time Transient Survey (CRTS) and SDSS Stripe 82. 
By comparing the expected Slepian Wavelet Variance from 
an AR(1) p rocess to that observ ed for real quasars in S82 
and CRTS, IGraham et al.l (I20l4 find that there is a clear 
indication that the AR(1) process fails to characterize AGN 
variability on short timescales. 

While one cannot ignore the success of the AR(1) pro¬ 
cess at modeling the time variability behavior of AGN on the 
timescales probed by MACHO, OGLE, SDSS, etc., we cau¬ 
tion that it may be impossible to observe deviations from the 
AR(1) process without probing the short timescales avail¬ 
able through Kepler. 


12 CONCLUSIONS 

Individual Kepler AGN light curves show a wide range 
of behavior that can be loosely grouped into 5 categories: 
stochastic-looking, somewhat stochastic-looking-|-weak os¬ 
cillatory features, oscillatory features dominant, flare fea¬ 
tures dominant, and not-variable. Some light curves appear 
to transition from one variability state to another, suggest¬ 
ing that AGN light curves may not be stationary in the sta¬ 
tistical sense (Sec. [ 6 ] and Fig. [3)). We estimate PSD slopes 
ranging from 7 = 1.7 to 7 = 3.1 with 3-8 objects having 
slopes close to that of the AR(1) process (7 = 2.0). Seven 
of the remaining objects have slopes significantly steeper 
than 7 = 2.5 suggesting strongly correlated non-AR(l) light 
curves, while two have 7 < 1.9. Four out of 5 of the AGN 
that exhibit oscillatory features in the light curve have PSD 
slope close to that of the AR(1) process, while the AGN that 
exhibit flares in the light curve are poorly fit by the DPL 
model (Sec. [To] and Table [3|). 

A broad dip feature can be observed in the structure 
function on timescales ranging from 10-100 d in the rest 
frame of 12 of the 20 AGN. This dip feature in the struc¬ 
ture function corresponds to increased correlation on these 
timescales in the light curve, i.e. flux measurements on 
timescales ranging from 10-100 d are closer than they 
should be. Due to the wide range of redshifts of these ob¬ 
jects, this timescale range is not identical in observed frame 
of these objects, suggesting that this feature is probably not 
caused by some sort of instrumentation issue (see Fig. HOI) . 
Some AGN exhibit varying structure function slopes, imply¬ 
ing that the value of 7 in the DPL model varies on different 
timescales (see Fig. 1131) . Lastly, not all AGN show variability 
at the levels probed by Kepler, confirming previous findings 
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that flux variability is seen in most but not all AGN (see 

Fig. [T2|. 

There is no clear relationship between the PSD slope 
and object type or redshift, although the sample size is 
too small to draw any deflnitive conclusion. We conclude 
that while enough of our light curves are consistent with 
an AR(1) process for the damped random walk model to 
be an appealing choice-especially for poorly sampled light 
curves-there are enough interesting features present in the 
light curves to warrant a more detailed analysis. 

The AGN light curves seen in Kepler suggest that AGN 
variability is a very complex phenomenon with individual 
light curves looking very different. We will perform a recali¬ 
bration of the Kepler AGN light curves to determine which 
light curve visual features are persistant. We will apply more 
sophisticated statistical models drawn from the fi eld of time 
series analysis such as the ARMA process model iHamilto^ 
(Il994h . The AR(1) process or DRW model can determine ‘if’ 
an AGN is varying, but is not helpful in determining ‘why’ 
the AGN is varying. 
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